An Evidence Based Search Method For Gravitational Waves From Neutron Star 

Ring-downs 
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The excitation of quadrupolar quasi-normal modes in a neutron star leads to the emission of 
a short, distinctive, burst of gravitational radiation in the form of a decaying sinusoid or 'ring- 
down'. We present a Bayesian analysis method which incorporates relevant prior information about 
the source and known instrumental artifacts to conduct a robust search for the gravitational wave 
emission associated with pulsar glitches and soft 7-ray repeater flares. Instrumental transients are 
modelled as sine-Gaussian and their evidence, or marginal likelihood, is compared with that of 
Gaussian white noise and ring-downs via the 'odds-ratio'. Tests using simulated data with a noise 
spectral density similar to the LIGO interferometer around 1 kHz yield 50% detection efficiency and 
1% false alarm probability for ring-down signals with signal-to-noise ratio p = 5.2. For a source at 
15 kpc this requires an energy of 1.3 x 10 _5 Mqc 2 to be emitted as gravitational waves. 

PACS numbers: 02.50.Cw, 04.80.Nn, 07.05. Kf, 95.55.Ym, 97.60.Jd 



I. INTRODUCTION 

A possible mechanism for the emission of gravita- 
tional waves from neutron stars is the excitation of non- 
radial quasi- normal modes (QNMs) This excitation 
could be caused by the disruption associated with pulsar 
glitches Q or from flaring activity in soft 7-ray repeaters 

The frequencies and damping times of the QNMs de- 
pend strongly on the neutron star equation of state 
(EOS) and for the more dominant /-modes these are 
thought to lie somewhere in the region of 1 to 4 kHz 
and 50 to 500 ms, respectively. Andersson & Kokkotas 
[H have shown how the mass and radius of a neutron star 
may be constrained by gravitational wave observations of 
the QNM frequencies and decay times. Conversely, if the 
EOS of the neutron star were known with some precision, 
it would be possible to compute the expected decay times 
and frequencies of the QNMs. This would provide a well 
constrained waveform model for aiding the identification 
of a potential gravitational wave signal following a neu- 
tron star ring-down. However, the behaviour of matter 
at the densities found in neutron stars is not well un- 
derstood and there exist many different models for the 
neutron star EOS. It is, therefore, necessary to develop 
techniques which are robust to the uncertainties in these 
models. 

While the gravitational wave emission from neutron 
star QNMs is expected to be weak (inducing typical 
strain amplitudes of ~ 10~ 24 ), their detection is further 
hampered by the presence of instrumental glitches that 
can closely resemble short-duration gravitational wave 



signals. 

In this work, we demonstrate the feasibility of applying 
Bayesian inference to the robust detection of neutron star 
gravitational wave ring-downs through a process of model 
selection. We highlight how the methodology may be 
extended to include a more realistic range of glitches and 
developed into a multi-detector search. 



II. PERTURBED NEUTRON STARS & 
ASTEROSEISMOLOGY 

When the solid crust of a neutron star is severely 
disrupted or cracked, some of the stored elastic en- 
ergy is channeled into the oscillatory modes of the star. 
Quadrupolar excitations will then be strongly damped 
by gravitational wave emission The different modes 
may be labelled according to the spherical harmonic in- 
dices I and m which describe the angular dependence and 
number of nodes. 

The fundamental fluid mode, or /-mode (as first shown 
by Kelvin for the case of a non-rotating, uniform density 
star) has angular frequency 
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where M and R are the stellar mass and radius, respec- 
tively. This is also a reasonable estimate for more realis- 
tic equations of state [f| and we see that, for non-radial 
modes with 1 = 2, the /-mode pulsations have frequen- 
cies ~ 2 kHz taking the fiducial values M = IAMq and 
R = 10 km. Other modes, such as the pressure (p) and 
space-time (w) modes, have considerably higher frequen- 
cies than this. Gravitational wave interferometers like 
GEO600 @, LIGO Q and VIRGO @ are more sensi- 
tive at lower frequencies making the /-mode the most 
favourable for a gravitational wave search. The ring- 
down timescale, r/, is given by the ratio of the oscillation 
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FIG. 1: The /-mode frequencies & damping times for the 
variety of different equations of state (indicated by symbol 
shape) and the extremal neutron star masses (indicated by 
symbol weight) considered in [Io| . 



energy to the total power emitted as gravitational waves 
Q. This yields 77 - R(R/M) 3 , where M and R are the 
mass and radius of the neutron star, respectively. 

So, given some EOS, which defines the ratio M/R, it 
is possible to calculate the exact damping times and fre- 
quencies of the QNMs. Andersson & Kokkotas [J] do 
precisely this for a variety of equations of state and es- 
tablish empirical relations for the /-mode damping time 
and frequency. On the other hand, we can consider the 
inverse problem and use gravitational wave observations 
of the /-mode frequency and damping time to constrain 
the neutron star EOS. Indeed, there has already been an 
attempt to use the electromagnctically observed frequen- 
cies of potential torsional mode oscillations in two SGRs 
to constrain their equations of state @ . 

Here, our interest is in using what little is known about 
the neutron star equation of state to inform us with re- 



gards to sensible gravitational wave waveforms to search 
for. Fig. Q] shows the results of calculations by Bcnhar 
et al [10] of the /-mode frequency and damping time, 
using several realistic EOSs and the extremal neutron 
star masses. M max refers to the maximum neutron star 
mass allowed by the EOS. The reader is directed to 
and the references therein for descriptions of the different 
equations of state. 



Examining Fig. [T] we see that the equations of state 
considered in yield typical /-mode frequencies of 
Vf ~ 1.5 — 3 kHz and damping timescales of 77 ~ 
50 — 400 ms. Later, we use these ranges to set sensible 
limits of our priors for each parameter. 

Gravitational wave emission 

Following Thorne [l| , we model the gravitational wave 
strain amplitude at the Earth from the n-th QNM as: 

r (t-*o) i 

ho sin [oj n (t — to)] e L T ™ J for t > to 
hil)={ (2) 

otherwise, 

where ho, u> n and t„ are the initial amplitude, angular 
frequency and characteristic damping time of the signal, 
and to is the start time of the signal. We take n = to 
represent the /-mode. 

Furthermore, ho, luq and r are related to the total grav- 
itational wave energy via the following relation Q , 



c 3 D 2 / i\2 
W = ~4G~ \ hoUjQT2 



(3) 



where D is the distance to the source. We can, there- 
fore, write down an expression for the the expected initial 
amplitude for /-mode ring-downs using fiducial parame- 
ter values: 



ho 



1.6 x 10-- ( E r ? ) (-^-) ~ 1/2 (-OL-) 
yiQ- 11 Mqc 2 J V200 ms/ V2 kHz/ 

r 



D 



15 kpc 



(4) 



Search triggers & QNM excitation 



There are several ways to generate stellar pulsations 
including close encounters with orbital companions and 
accretion (e.g., of comets), the ringing of proto-neutron 
stars following supernovae, soft 7-ray repeater (SGR) 
flaring activity and crustal disruption due to pulsar 
glitches. 



Our aim is to search for gravitational waves from the 
latter two mechanisms. This is because both pulsar 
glitches and the SGRs provide an observable electromag- 
netic counterpart (which may be used to 'trigger' the 
search) and occur with a frequency making them practi- 
cal for a triggered search. The most prolifically glitching 
pulsar, PSR J0537-6910, has a glitch rate of ~ 4 year -1 
[llj |. whilst other regular glitchers like the Crab pulsar 
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and Vela have rates of order ~ 0.5 year -1 [T3.[T3T|. SGRs 
are seen to emit bursts at rates of a few to 10 year -1 
[T3 |. whereas hyperflares are much rarer events occuring 
maybe once a decade. Other potential candidates seen 
to glitch are the anomalous X-ray pulsars (AXPs). 



Pulsar glitches 

Pulsar glitches are observed as sudden irregularities 
in the rotation rate of pulsars and are characterised by 
a step increase in rotation frequency. The characteris- 
tics of glitches vary between pulsars. For example, some 
show step changes in the spin-down rate and others ex- 
ponential recoveries to pre-glitch parameters. The exact 
characteristics give clues to the underlying mechanisms 
causing the glitch. 

The mechanisms responsible for pulsar glitches are still 
unclear, but there are two main candidates to explain the 
underlying process and some of the differences between 
glitches. For older pulsars it seems glitches are likely 
caused by a dramatic decoupling between the star's solid 
crust and superfluid interior [15fl . Glitches in the young 
Crab pulsar are thought to be associated with a reconfig- 
uration of the solid crust [l6j : the spin-down reduces the 
centrifugal force and the crust reaches breaking strain. 
The ensuing relaxation of ellipticity will cause a sudden 
change in the moment of inertia, producing the observed 
glitch, and the crustal rupture will set up a starquake, 
hopefully causing /-mode excitations. The reality of the 
mechanism is likely to be very complex and may be a 
combination of the two. 

For the two different glitch models the amount of en- 
ergy released can be estimated in different ways as shown 
in van Riper et al. (1991) [l7|. They assume all the en- 
ergy released goes into heating the star, whereas we will 
make the assumption that it goes into exciting quasi- 
normal modes. For the angular momentum exchange 
model (thought to be the most probable explanation for 
the Vela pulsar glitches) the amount of energy released 
depends on the angular momentum exchanged between 
the superfluid interior and crust AJ ~ I Aft, where Af2 
is the angular frequency change from the glitch. The en- 
ergy change is then AE — AJf2i ag , where fii ag is the 
lag frequency between the superfluid and crust, with an 
estimated range of values of 1-100 rads -1 (or possibly 
< O.lrads -1 ) [17]. For the largest Vela pulsar glitch, 
with a fractional frequency change of Ail/57 = 3.1xl0 -6 
[TH . this gives a A J ~ 2 x 10 34 J giving an energy re- 
lease of AE ~ 10" 13 - 10 -11 M Q c 2 for the range ^i ag w 
1-100 rads" 1 . 

For starquake driven glitches the energy released is 
given in Rcf. 17] as AE k /xT4 rus te max eq Uake , where 
£quakc = Afl /ft is equivalent to the relative change in mo- 
ment of inertia, /i is the mean shear modulus of the star, 
V crus t is the volume of the crust (where {J.V crus t ~ 10 41 J), 
and e max is the maximum deformation from equilib- 
rium the crust can withstand without breaking (given 



in Ref. [I?) as e max < 10 -2 although this could vary 
somewhat). Assuming this e max and taking the largest 
Crab pulsar glitch, where Afl/fl ~ 8 x 10 -8 [l]|, we 
get an energy release of AE ~ 5x 10 -16 M Q c 2 . If the 
starquake mechanism can provide similar fractional fre- 
quency changes to a neutron star to those seen in the 
Vela pulsar during glitches, then this mechanism could 
still be a valuable potential source. 



Soft y-ray repeater flares 

Soft 7-ray repeaters are high energy transient sources 
with typical photon energies of 10 — 30 keV and similar 
burst characteristics from one event to the next, although 
they are also seen as quiescent X-ray sources. These 
objects are identified as highly magnetised (B w 10 14 
Gauss) neutron stars, or 'magnetars'. They are occasion- 
ally seen to emit giant flares, or hyper-flares, which have 
thousands of times the luminosity, and harder spectra, 
than the regular bursts. The hyperflares are thought to 
occur when magnetic field becomes twisted and causes 
a catastrophic reconnection, inducing tectonic activity. 
The field annihilation/reconnection in seismic faults is 
responsible for the observed 7-ray emission and, again, 
we expect the crustal disruption to excite the QNMs of 
the magnetar [20| . 



III. BAYESIAN MODEL SELECTION 

Before moving on to describe the search in detail, we 
outline some ideas behind Bayesian model selection. 

Bayes' theorem describes how to assign a posterior 
probability to some parameterised model Mj, given a set 
of data or observations D and some background informa- 
tion / which determines the hypothesis space {Mi}: 

p( M, |a/) = SM. (5) 

Here p(D\Mi,I) is the marginal likelihood or evidence 
and represents the influence of the data on our belief in 
Mi 1 ; p(M l \I) is the prior probability of model Mi and 
describes our state of belief in Mi, preceding examina- 
tion of the data D; p(D\I) ensures that the posterior is 
correctly normalised. 

Our desire here is to detect a gravitational wave of 
a known shape (see Eq. ^ but with unknown parame- 
ters within some range in noisy data. The obvious ques- 
tion we might ask is, 'does this data contain a ring-down 



The name marginal likelihood reflects the fact that we have 
marginalised over the parameter values that are associated with 
Mi. 
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gravitational wave?'. We denote the data by D and the 
gravitational waveform model by Mqw- 

The posterior probability p(Mi\D,I) then tells us the 
degree of belief to assign to the model Afj. To properly 
normalise the posterior, however, we must marginalise 
over the entire hypothesis space: 



p(d\i) = J2p(m 1 \Dp(d\m 1 ,i). 



(6) 



If the system is sufficiently well understood, it is possi- 
ble to enumerate all possible models Mi and the posterior 
probability of any one model can be calculated directly 
from Bayes theorem. However, this cannot always be 
done and it often makes more sense to evaluate the prob- 
ability of one model relative to another. Such compar- 
isons are performed via the odds ratio: 



12 = 



p(M 2 \D,I)- 



(7) 



For a gravitational wave search we might choose Mi to 
be the proposition that the noisy data contains a gravita- 
tional wave, and M 2 to be the proposition that the data 
only contains detector noise. 

Substituting the right hand side of Bayes' theorem for 
the posteriors in Eq. we see that the normalisation 
term p(D\I) drops out and we are left with 



pjM^I) pjD^hJ) 
p(M 2 \I)p(D\M 2> iy 



(8) 



The first term, the prior odds, is the ratio of the prior 
probabilities for each model. Typically, we assume com- 
plete naivete and set the prior odds equal to unity. The 
second term is the ratio of the evidences from each model 
and is called the Bayes factor. Clearly a large value of 
the Bayes factor indicates a strong preference for M\. 

The evidence is computed by integrating the likelihood 
p(D\6, Mi, I) over all model parameters 6 and weighting 
by the prior on those parameters, p(8\I), leading to the 
alternative name 'marginal likelihood': 



P (D\Mi,i) = [ p(e\Mi,i)p(D\e,Mi,i)de. 

Je 



(9) 



So given some competing models M\ and M 2 , we can 
evaluate the evidences p(D\Mi, I), p(D\M 2 ,I) and as- 
suming prior odds of unity, use the odds ratio to decide 
which is most likely, given the data D. 



A. Application & choice of models 

Here, 'model' shall refer to a class of descriptions for 
the data. An example is 'the data contains a ring-down 



gravitational wave signal in addition to noise'. Note that 
we have defined the generic shape of the data (a noisy 
ring-down) but not any parameter values. A particular 
signal with a specified set of parameter values is called a 
'template'. In this way, a model defines a set of templates 
with parameter values determined by the priors in that 
model. Additionally, when we talk about the 'evidence 
for the model', we are referring to the marginal likelihood 
for that model, i.e., p(D\M, I). The total evidence of the 
hypothesis space p{D\I) is eliminated through the use of 
the odds ratio. 



Mi: Ring-down waveform & Gaussian white noise 



The expression for the ring-down waveform h(t) is 
given by Eq. O We assume that the noise n(t) is white 
and Gaussian over a sufficiently broad band and that the 
data stream is given by 



d{t) = h(t) +n(t), 



(10) 



where the noise n(t) has zero mean and variance a\. 
However, to simplify data conditioning, we work with 
the power spectral density, 



D(w) = \h(uj)\ 2 + \n(uj)\ 2 + 2\h(uj)n(uj)\, 
where h(u>) is the Fourier transform of h(t), 



\h(u)\ 



h T 



^/l + ^-wo) 2 



(11) 



(12) 



so that our parameter space is given by 6 = {ho, too, r}. 
Notice that working with the power spectral density has 
the effect of pre-marginalising over the start time of the 
signal, to with a uniform prior. The power spectral den- 
sity is estimated from fast Fourier transforms (FFTs) of 
consecutive segments of the time series. This yields a 
spectrogram with time bins indexed by i and frequency 
bins indexed by j. 

The likelihood of obtaining power Dij in the i-th time 
bin at the j-th frequency (i.e., the i,j-th pixel), given 

is a non- 



a template with signal power Si 



central x 2 distribution with two degrees of freedom and 
non-centrality parameter equal to the power from the 
template, i.e. 
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exp 



where afj is the variance of the Fourier components in 
that pixel and Iq is the zeroth order modified Bessel 
function of the first kind. If the data is Gaussian and 
white, the power spectral density can be normalised such 
that the Fourier components follow Gaussian distribu- 
tions with zero mean and unity variance (i.e., cr|- = 1). 

To calculate a single odds ratio in each time bin, we 
require the joint probability across frequencies, 

p({D}\e,Mi,I) = Y[p(D j \0,M 1 ,I), (14) 

3 

where we have dropped the time bin index i for notational 
convenience. Finally, we adopt independent, uniform pri- 
ors on the parameters 6 = {}iq,luo,t} with ranges on loq 
and r based loosely on the expected values in the litera- 
ture. For the amplitude ho, however, there is little to be 
gained by restricting the prior range and the prior on ho 
is taken to run from zero to some arbitrarily high value. 
The joint prior is then 

p(e\M u I)=p(h \M 1 ,I)p(co \M 1 ,I)p(T\M 1 ,I) (15) 

where 

{g _g . for 9 min < 9 < # max 
(16) 
otherwise, 

and 9 is any one of the parameters in model 1. 



M%: Gaussian white noise only 

M2 is our null detection hypothesis: the data is mod- 
elled as white Gaussian noise, without any gravitational 
wave signal, so that d{t) = n(t). Again, we work in terms 
of the power in each spectrogram pixel D,j . If the spec- 
trogram has been normalised such that the individual 
Fourier components are normally distributed with mean 
zero and variance of unity, we know that the power in 
the J-th frequency bin follows a central % 2 distribution 
with two degrees of freedom, 

p(A,|M 2 ,/) = ^expj-0 j, (17) 

where tr| is the variance of the Fourier components in 
the j-th frequency bin. If it is possible to estimate the 
variance in each frequency bin (from an off-source piece 
of data, for example), then the result is also valid for 



1 

coloured noise. 

Notice that we have arrived at the evidence from Mi 
without making any mention of parameterisation, priors 
or marginalisation. This can also be derived in a purely 
Bayesian context by considering the likelihood of a pixel 
power Dij , given a template power Sij (Eq. I13[) . In the 
case of a model where there is no contribution to the 
power from a gravitational wave, we know a priori that 
Sij — 0. That is to say 

p{S ij \M 2 ,I) = 6(S ij ), (18) 

where S is the Dirac delta function. If we now marginalise 
the likelihood given by Eq. 1 131 over power using this prior, 
we arrive at Eq. [T7l 



Thresholding O12 

O12 is the ratio of the posterior probabilities for each 
model, so it might seem sensible to choose O12 > 1 to 
indicate a preference for Mi over M^. While this is true, 
such a threshold for gravitational wave detection neglects 
the role of our prior odds and the need for an acceptable 
false alarm rate. 

By setting the prior odds equal to unity, we are saying 
we believe a priori that both models are equally proba- 
ble. Even if we truly were that ignorant, the influence of 
the data through the Bayes factor may cause the odds to 
fluctuate around some mean value away from unity. In- 
stead, we search for excesses from the mean 'off-source' 
(zero-signal) odds to indicate a preference for our gravi- 
tational wave model. Alternatively, it would be straight- 
forward to estimate the prior odds using an off-source 
sample. The value of the prior odds could then be cho- 
sen such that an odds ratio of unity corresponds to a false 
detection probability of 0.5. Ultimately, the odds thresh- 
old, denoted Othresh, can be set according to the results 
of large numbers of trials and so the overall normalisation 
is relatively unimportant. 

Finally, the magnitudes of the variations of the odds 
ratios mean that it is most natural and convenient to 
work with the base 10 logarithm of the odds ratio. For 
clarity, log 10 O12 refers to the base 10 logarithm of the 
odds ratio between models Mi and M^. 
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B. Algorithm &: example 

We now consider an example using data synthesised 
in Matlab 2 to illustrate the above principles. First, the 
outline of the algorithm is as follows: 

1. Estimate the variance of the noise using some 
stretch of data away from the time of an expected 
gravitational wave signal (i.e., off-source). 

2. Construct the power spectral density of discrete 
time segments of data centred around the expected 
gravitational wave signal (i.e., on-source) to create 
a time- frequency map of power (spectrogram). 

3. Normalise the power spectral density so that, in the 
absence of a signal in the data, the power in a given 
frequency bin follows a central % 2 distribution with 
two degrees of freedom. 

4. Compute the evidences p(D\M\, I), p(D|M 2 , /) for 
each model in each spectrogram time bin. 

5. Assuming prior odds of unity, evaluate the odds 
ratio O12 in each time bin. An excess in the odds 
ratio indicates a preference for Mi and, therefore, 
a potential detection. 

To characterise the signals used to test the algorithm, we 
define the signal-to-noise ratio as 



S{v) 



(19) 



where S(y) is the one-sided noise spectral density. In the 
case of Gaussian white noise, S(u) is given by 
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FIG. 2: The normalised power spectral densities of the ring- 
down (dashed line) and sine-Gaussian (solid line) injections. 



Injection 


Parameter 




Value 


RDI 


Injection time 


to 


10s 




Initial amplitude 


h 


7.0 x 10~ 21 




Central frequency 




2 000 Hz 




Decay time 


T 


0.21s 




SNR 


P 


20.4 


SGI 


Injection time 


to 


80s 




Initial amplitude 


h 


3.72 x 10~ 21 




Central frequency 


vo 


2 000 Hz 




Decay time 


T 


0.21s 




SNR 


P 


20.4 



Is 



(20) 



where ct 2 is the variance of the time series data and f s 
is the sampling frequency. It is also useful to define the 
root-sum-squared amplitude 



+oo \ 1/2 

\h(t)\ 2 dt 



(21) 



To demonstrate the operation of the algorithm, we in- 
ject a ring-down signal into 100 s of Gaussian white noise 
with amplitude spectral density 10~ 22 Hz -1 / 2 . To inves- 
tigate the response to a typical instrumental glitch which 
closely resembles our target waveform, we also inject a 
sine-Gaussian signal of the form 



h{t) = h sin [oj (t - t )} e~ (t " to)2/r2 



(22) 



1 http : //www.mathworks . com/products/mat lab 



TABLE I: Injected signal parameter values. 



where, hq is the maximum amplitude of the signal, uj is 
the angular frequency and r is the decay time. Table U 
shows the parameter values used to generate the injec- 
tions and Fig. [2] shows the power spectral density of each 
signal, calculated from the noisy data and normalised so 
that the noise follows a central \ 2 distribution. It is the 
job of the algorithm to detect and differentiate between 
these two signals, only producing a candidate detection 
when the ring-down is present. The prior ranges used for 
each parameter are shown in Table HT1 



Fig. [3] shows the odds ratio in each spectrogram time 
bin. The ring-down injection at t = 10 s is strongly 
detected. Notice, however, that the sine-Gaussian we 
have injected to mimic an unwanted instrumental glitch 
is also detected, albeit with an odds ratio lower than 
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TABLE II: Parameter prior ranges. The prior distributions 
are taken as uniform over these ranges. 
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FIG. 3: The log-odds in favour of a ring-down plus noise, 
versus the possibility that the power is due to noise alone. 
The first peak corresponds to the ring-down injection, the 
second to the sine-Gaussian injection. 



that for the ring-down injection. To address this issue, 
we make a straightforward extension to the odds ratio to 
consider multiple hypotheses. 



C. Multiple hypothesis extension 

To handle the possibility of sine-Gaussian instrumental 
glitches, we rewrite the posterior in the denominator as 
the sum of the probability of the noise model M 2 and a 
model for sine-Gaussian glitches, M%. This gives 



p(M_|£>, I) = p(M 2 \D, I) + p(M 3 \D, I) 



(23) 



where M_ denotes the proposition that the data does not 
contain a ring-down gravitational wave and the ellipsis is 
to emphasise the fact that this null-detection hypothe- 
sis M_ may be further extended to include additional 
models for instrumental glitches. Similarly, 



p(M+\D,I)=p(M 1 \D,I) 



(24) 



where, again, it is straightforward to include additional 
signal models if desired. The result is that we are left 
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FIG. 4: The log-odds in favour of a ring-down plus noise, 
versus the possibility that the power is due to noise alone or 
a sine-Gaussian glitch plus noise. The second (sine-Gaussian) 
peak seen in Fig. [3] is now seen as a dip. 



with a comparison of the probability in favour of a grav- 
itational wave with the probability of an instrumental 
glitch or that of the noise model, thus making maximal 
use of any knowledge we might have regarding transient 
features in the data. 

Using M_ now as the alternative hypothesis, we obtain 
a new odds ratio O123: 



O123 — 



p{M x \D,I) 



p(M 2 \D,I)+p(M 3 \D,iy 



(25) 



Assuming the prior on each model is identical, the pos- 
terior probabilities in O123 may be replaced by the evi- 
dences p(D\Mi, I) for i = 1, 2, 3. 

Fig. 0] shows the result of using this new expression for 
the odds ratio and the same data from the previous ex- 
ample. Since the parameterisation for the sine-Gaussian 
is identical to that of the ring-down, the priors in table 
HJare also used to evaluate p(Ms\D, I). There is no fun- 
damental reason for using the same prior ranges for both 
models. Indeed, a realistic application would most likely 
have very different priors for the parameters in different 
models even if the parameterisation was the same. Here, 
the priors in Table [TT| are used for both Mi and M3 for 
computational simplicity. 



The incorporation of the alternative model eliminates 
the previous problem of detecting sine-Gaussian signals 
with high odds ratios and, in fact, the odds ratio now 
strongly prefers the sine-Gaussian model to the ring- 
down model in the presence of the sine-Gaussian. Also 
note that the size of the peak indicating the ring-down 
has diminished quite substantially but is still clearly vis- 
ible above the background. This reduction is due to a 
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non-zero contribution from p(M^\D, I) in the denomina- 
tor of Oi23- Unless they are mutually exclusive in some 
way, the inclusion of additional models in the denomina- 
tor of the odds ratio will generally increase the robustness 
of the search, at the cost of sensitivity. 

IV. PERFORMANCE 

We now investigate the performance of the algorithm 
by considering both formulations of the odds ratio, com- 
paring the relative merits of each. 

A 'false alarm' is defined as any unwanted transient 
event which causes the odds ratio to cross the thresh- 
old and generate a candidate detection. For white noise, 
false alarms are caused by spikes in the noise amplitude. 
The false alarm probability is calculated from the frac- 
tion of time bins in a large sample for which the odds 
ratio crosses the threshold when there is no ring-down 
signal present. 

We investigate the response to 'off-source' data (i.e., 
white noise with no injections) by combining the results 
of ten 500 s spectrograms with a 1 s time resolution to 
give 5 000 off-source time bins. To evaluate the sensitivity 
of the search, ring-down signals of a constant signal-to- 
noise ratio are injected every other second into 500 sec- 
ond segments of white noise, synthesised in Matlab. This 
is then used to construct a 500 s spectrogram with 1 s 
time-resolution and 0.5 Hz frequency resolution for each 
set of injections, leading to 250 signal injections for each 
value of the signal-to-noise ratio. The fact that there are 
spectrogram time bins with no signal injection helps to 
prevent contamination from adjacent bins. We vary the 
signal-to-noise ratio through the value of ho only and al- 
ways compare signals of equal bandwidth and at the same 
frequency, with the frequency held constant at 2 kHz and 
the decay time at 207.5 ms. The noise amplitude spec- 
tral density is 10 _22 Hz -1 / 2 , representative of the LIGO 
noise floor at these frequencies. 

When we make use of the glitch catalogue in the ex- 
panded odds ratio O123 the objective is to be robust 
against unwanted glitches as well as spurious noise ef- 
fects. Therefore, we also define a false alarm probabil- 
ity due to sine-Gaussians. This is found in exactly the 
same way as the sensitivity to ring-downs but using sine- 
Gaussian injections. The false alarm probability due to 
sine-Gaussians, therefore, is simply the fraction of the 
250 injected sine-Gaussians which cause the odds ratio 
to exceed the threshold. 



A. O12: ring-down vs noise 

We take our desired false alarm probability here to 
be 1%. Fig. [5] shows the false alarm probability as a 
function of the log odds ratio threshold, log 10 Othresh, 
and we see that a false alarm probability of 1% is given 

by log 10 Othresh = 4. 




-2-101234 



FIG. 5: The false alarm probabilities in Gaussian white noise 
for different log 10 O12 thresholds. Error bars represent la 
Poissonian standard errors. The solid vertical line marks the 
threshold applied when using log 10 O12. 



Fig. [5] shows receiver operating characteristic (ROC) 
curves for Oyi- ROC curves are plots of sensitivity (the 
probability of detecting what we are looking for) as a 
function of false alarm probability (the probability of 
claiming a detection due to a spurious noise event or in- 
strumental glitch) for a given strength signal. The differ- 
ent false alarm probabilities correspond to different val- 
ues of the detection threshold, Othresh, and are found by 
reading the appropriate values from Fig. [5] When there 
are no injected signals, the only events to cause the odds 
ratio to cross the threshold are false alarms and the false 
alarm probability should be equal to the detection prob- 
ability. 



Although the ROC curves give some indication of how 
the sensitivity varies with injected signal strength, it is 
also helpful to examine efficiency curves, where the frac- 
tion of detected ring-downs is plotted against the injected 
signal strength, for a given threshold. Fig. [7] shows the 
detection efficiency using log 10 Othresh = 4 with corre- 
sponding false alarm probability of 1%. 



The algorithm's performance is best summarised by 
the signal-to-noise ratio required to give a detection prob- 
ability of 50% while maintaining the desired false alarm 
probability. A signal-to-noise ratio p — 5.2 is required 
for 50% ring-down detection, corresponding to an initial 
ring-down amplitude at Earth of /to = 1.8 x 10~ 21 at 
current LIGO sensitivities. 

Equation [3] on page relates ho to the distance to the 
source and the energy emitted as gravitational waves. 
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FIG. 6: Receiver operating characteristic curves for On. 
False alarm probabilities are those from amplitude fluctua- 
tions in Gaussian white noise. 
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FIG. 7: The detection efficiency obtained for log 10 O12 using 
a threshofd log 10 Othrcsh — 4. Error bars show la standard 
Possonian errors. 

Using this equation and assuming the fiducial distance 
of 15 kpc, the energy required to generate this initial 
ring-down amplitude at the Earth is 1.3 x 10~ 5 M@c 2 . 
Conversely, when we assume the fiducial energy in grav- 
itational waves of 10 _11 Mq, the distance to the source 
must be 13.2 pc for this amplitude. These results are 
summarised in Table Hill 



B. O123: Ring-down vs noise or sine-Gaussian 
glitch 

We now evaluate the performance of the algorithm 
using the more robust comparison between ring-downs, 



FIG. 8: The mean value of log 10 O123 for varying strengths 
of signal-to-noise ratio. Results of ring-down injections are 
shown by the dashed curve, while the solid curve shows the 
results of the sine-Gaussian injections. Error bars represent 
la standard errors in the mean values of log 10 O123. 

white noise and sine-Gaussians. 

The response of O123 to sine-Gaussians is compared 
with the response to ring-down injections in Fig. [H] The 
horizontal axis shows the signal-to-noisc ratio of each 
type of injected signal (ring-down or sine- Gaussian) and 
the vertical axis gives the mean value of log 10 O123, cal- 
culated from 250 injections of each signal type. The re- 
sponse can be separated into three regions of signal-to- 
noise ratio, p: 



p < 3 : there is no response to sine-Gaussians or ring- 
downs and any signal present is indistinguishable 
from noise. 

3 ^ p 5 : the algorithm begins to respond to the ring- 
down injections. The sine-Gaussian injections gen- 
erate a very slight rise in the log odds. 

p ^ 5 : the odds in favour of a ring-down begin to grow 
for the ring-down injections and rapidly falls off for 
the sine-Gaussian injections. 

We use the intermediate region (i.e., where 3 < p < 5) 
to assume a worst-case scenario where all sine-Gaussian 
glitches have a signal-to-noise ratio p ~ 4 and compute 
the false alarm probability from sine-Gaussians for dif- 
ferent odds ratio thresholds. The results are shown in 
Fig. [3 



We now find that the threshold required to give a 
false alarm probability of 1% in Gaussian white noise 
is log 10 Othresh = 0.84. This corresponds to a false 
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FIG. 9: Top panel: false alarm probabilities and correspond- 
ing log 10 O123 thresholds for Gaussian white noise. Bottom 
panel: false alarm probabilities and corresponding log 10 O123 
thresholds for sine-Gaussians with signal-to-noise ratio p ~ 4. 
In each case, the solid vertical line indicates the threshold 
required for a 1% false alarm probability in Gaussian white 
noise. Error bars show la Poissonian standard errors. 



alarm probability of 10% when the data contains a sine- 
Gaussian. We note two points here: 

a) we have assumed the presence of sine-Gaussian 
glitches a priori. While it is appropriate to as- 
sume the presence of Gaussian white noise, we do 
not have a population model for the sine-Gaussian 
glitches. A more informative estimate of the proba- 
bility of mistakenly detecting a sine-Gaussian glitch 
should fold in the effects of such a population model 
through the prior on M3. In this work, we are 
more interested in demonstrating the inherent abil- 
ity of the algorithm to discriminate between similar 
waveforms and this will not be considered. 

b) Only sine-Gaussians with signal-to-noise ratios of 
p ~ 4 are considered. Away from this value, the 
odds ratio drops off rapidly so that the 10% false 
alarm probability can be regarded as an upper limit 
for sine-Gaussian glitches. 
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FIG. 10: Receiver operating characteristic curves for O123. 
Top panel: ROC curve with false alarm probabilities due to 
spurious noise events. Bottom panel: ROC curve with false 
alarm probabilities due to sine-Gaussian glitches. 



Because we have two possible sources of false alarms, 
sine-Gaussians and noise events, we produce ROC figures 
for each, shown in each of the panels of Fig. ITUl 

For comparison with the sensitivity of O12 Fig. 1111 
shows the detection efficiency for O123. We find that the 
signal-to-noise ratio required for 50% detection efficiency 
and 1% false alarm probability in Gaussian white noise 
is p 5Q% = 8.0. 



Again, it is important to consider these results in an 
astrophysical context. With the noise amplitude spec- 
tral density used for these investigations (10~ 22 Hz -1 / 2 ), 
a ring-down signal-to-noise ratio p = 8.0 corresponds to 
an initial amplitude ho — 2.9 x 10~ 21 Hz -1 / 2 . The grav- 
itational wave energy required to produce a ring-down 
signal with 50% detection probability at the Earth with 
this amplitude is 3.3 x 10 -5 Mqc 2 for a source at 15 kpc. 
Similarly, if we assume that 10 -11 Mqc 2 is emitted as 
gravitational waves , the source must lie at 8.2 pc. These 
results are summarised and compared with those from 
X2 in Table EU 
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Parameter 


Symbol 


Value Using Oi 2 


Value Using 


0123 


Log odds threshold 


log 10 Othresh 


4 


0.84 




False alarm probability (due to noise) 


p(-\N) 


1% 


1% 




False alarm probability (due to sine-Gaussians) 


p(-\SG) 


— 


10% 




Signal-to-noise 


p50% 


5.2 


8.0 




1111 lldl tall dill dill Ull L LLU.C 


u 50% 

"o 


i a v 1 n — 21 

± .o X 1U 


2.9 x 10" 


21 


GW Energy 


P 50% 
fi GW 


1.3 x 1(T 5 M0C 2 


3.3 x lO- 5 M c 2 


Range 


n 50% 
"GW 


13.2 pc 


8.2 pc 




Root-sum-squared strain 


1,50% 
rss 


4.1 x 10- 22 


6.6 x 10 - 


22 



TABLE III: Simulated sensitivity estimates for 50% detection efficiency and corresponding false alarm probabilities. p( — \N) is 
the probability of a false alarm from white noise for the given threshold log 10 Othresh- Similarly, p( — \SG) is the probability of 
a false alarm given a sine-Gaussian event. All estimates assume a ring-down frequency v = 2 kHz, decay timescale r = 207.5 
ms and a noise amplitude spectral density 10~ 22 Hz~ 1//2 . Eq^j IS calculated for a distance D — 15kpc and Dq^j 1S calculated 
for an energy in gravitational waves of 10 _11 A/qc 2 . 
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SNR of software injected waveform 

FIG. 11: The detection efficiency obtained for log 10 O123 us- 
ing a threshold log 10 Othresh — 0.84. The solid horizontal 
line indicates 50% detection efficiency. Error bars show lex 
Poissonian standard errors. 



holes [22| • The software can be used to compare the false 
alarm probability of the matched filtering method with 
our method by running it on simulated white noise for the 
range of parameters given in 

Table OH 4 . The template 
bank for the matched filtering was produced to give a 
maximum mismatch between adjacent templates of 1%. 

We define the false alarm probability to be the fraction 
of detection candidates generated by the algorithm whose 
signal-to-noise ratio crosses a given threshold. Fig. 12 
shows the false alarm probability as a function of signal- 
to-noise ratio threshold: for a false alarm probabilty of 
1%, we require a signal-to-noise ratio threshold of 5.85 
(compared with p = 5.2 for the evidence based search). 

It may seem surprising that the evidence-based search 
marginally out-performs matched filtering. It is impor- 
tant to remember, however, that the approaches are by 
no means equivalent: matched filtering searches for con- 
sistency with a given waveform, whereas the evidence- 
based approach searches for both consistency with the 
same waveform but also for inconsistency with a noise 
model. Since the noise model is highly sensitive to ex- 
cess power (which exponentially decreases its likelihood) , 
it does not seem unreasonable that there is a slight dis- 
parity between the approaches. 



C. Comparison to matched filtering 

The LIGO Algorithm Library (LAL) 21] and LALapps 
software repositories used for much of the gravitational 
wave data analysis within the LIGO Scientific Collabora- 
tion (LSC) 3 contain software for performing a matched 
filter based search for ring-down signals. Currently this is 
being used to search for ring-downs from perturbed black 



V. FURTHER EXTENSIONS 

We now highlight some outstanding issues and possible 
important extensions to this work which may be required 
to make the transition from 'proof-of-concept' to a use- 
ful search tool. Particularly, we address the issues of 



3 http://www.ligo.org 



4 The matched filtering code defines the ring-down in terms of 
frequency and quality factor, Q, where Q = nrf. 
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FIG. 12: The cumulative percentage of triggers produced by a 
matched filter based ring-down signal search code over a range 
of signal-to-noise ratios for data consisting of simulated Gaus- 
sian white noise. Error bars indicate la Poissonian standard 
errors. The solid vertical line indicates the p = 5.85 threshold 
for a 1% false alarm probability. 



constructing the catalogue of alternative hypotheses, the 
extension to a multi-detector analysis and the outline for 
a full analysis pipeline. 



Glitch classification 



in this work, we aim to compute the posterior probabil- 
ity for some gravitational wave model, given the interfer- 
ometer data and some background information. In the 
multi-detector case we still seek the posterior probability 
for some model but now using the information contained 
in each detector's output. Suppose then that we have 
some gravitational wave model Mqw and detector out- 
puts D\ and D 2 . We can immediately write down the 
posterior for Mqw using Bayes' theorem 



p(M GW \D 1 ,D 2 ,I) 



p{M GW \I)p{D u D 2 \M GW J) 
p(D 1 ,D 2 \I) 



where the gravitational wave model Mqw factors in the 
appropriate detector response functions for the source 
sky position and for detectors with uncorrelated output, 
the joint probabilities are the products of the individual 
probabilities. For N detectors 



p(M, 



GW 



|P},/)=p(M GW |J)n P(D ;{^' J) (27) 



and we can again eliminate the denominator by compar- 
ing the relative probabilities of different models via the 
odds ratio. Notice, however, that if the probability from 
one detector is very large but low in the other detectors 
it is still possible to get a high posterior across all the 
detectors. It would, therefore, be sensible to apply a cut 
to the data such that the odds must cross some thresh- 
old for each individual detector before considering the 
multi-detector case. 



It is likely that the inclusion of extra information re- 
garding glitches will increase robustness against instru- 
mental transients which would otherwise generate a can- 
didate detection event. For this to be effective, we re- 
quire some classification scheme for these instrumental 
glitches which would allow the construction of a cata- 
logue of glitch models as well as their associated priors. 

Essentially, what is required is an automated pattern 
recognition tool which could be 'trained' on off-source 
interferometer data and used to generate a generic li- 
brary of transient features. Fortunately, on-going detec- 
tor characterisation work aims to perform such an anal- 
ysis [H|. 

Multiple detector case 

A potential issue with this search is in fact the preva- 
lence of instrumental ring-downs already present in the 
interferometer data. Unless such features can be vetoed 
using known instrumental couplings, for example [24| . 
the only way to distinguish these from the targeted grav- 
itational wave ring-down is to search for coincidences be- 
tween multiple detectors. 

Here, we benefit again from the simplicity of the 
Bayesian formalism. In the single detector case outlined 



Analysis pipeline 

Finally, we present the outline for a planned analysis 
pipeline in Fig. 1131 Here, we assume a network of two 
detectors and that there is no correlation in the output 
from each. 



Upon reception of an external trigger, such as a pulsar 
glitch or soft 7-ray repeater flare, we would retrieve inter- 
ferometer data from near the event (on-source data). The 
length of data required will typically depend on the time 
resolution of the external trigger. Additionally, we re- 
quire data preceding the external trigger to provide some 
estimate of the background odds ratio and hence set the 
threshold for the desired false alarm probability. 

Following the acquisition of both on and off-source 
data for each detector, we transform to the Fourier do- 
main and construct spectrograms of each. From here, 
the methods outlined in this work can be used to com- 
pute the evidences for each gravitational wave and glitch 
model and the odds ratio can be computed for each de- 
tector using on and off-source data. 

If the odds in both detectors cross some appropriate 
threshold, then we can go on to calculate the multi- 
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FIG. 13: Planned multi-detector, evidence-based analysis 
pipeline. 



detector odds ratio. If this multi-detector odds ratio also 
crosses a threshold set from off-source data then we have 
a candidate detection. 

Alternatively, if the odds ratio from just one of the 
detectors in the network crosses the respective threshold, 
then there is probably good evidence for a glitch. Further 
examination of the odds ratios for each glitch model may 
reveal its nature and the information could conceivably 
be used to update the prior on that particular model in 
future studies. 



VI. CONCLUDING REMARKS 

We have examined the feasibility of using a Bayesian 
evidence discriminator in an externally-triggered search 
for gravitational waves produced by neutron star quasi- 
normal mode ring-downs. The evidence may be thought 
of as the total probability of the data given some model. 
By comparing the evidences for competing models via the 
'odds ratio', we can naturally make maximal use of prior 
information regarding the source as well as any informa- 
tion available from detector characterisation studies. 

This is particularly important in the single-detector 
case. Here there may be instrumental transients or 
'glitches' which would otherwise generate a detection 
candidate in algorithms which only look for consistency 
with the target model. We begin with a simple formu- 
lation of the odds ratio, O12, where we compare the 
probability of a noisy ring-down model with the prob- 
ability of noise alone. We find this yields a 50% de- 
tection efficiency for ring-downs with frequency 2 kHz, 
decay time 0.21s and signal-to-noise ratio p = 5.2 for 



a 1% false alarm probability in Gaussian white noise of 
amplitude spectral density 10 -22 Hz -1 / 2 . To obtain the 
same false alarm probability using matched filtering we 
require a signal-to-noise ratio threshold p — 5.85. As de- 
scribed in section HV CI this disparity is most likely due to 
the fundamental difference between using the odds ratio, 
which compares the evidences from two or more compet- 
ing models and matched filtering, which only searches for 
consistency with a single target waveform. 

Assuming 10 _11 Mqc 2 is emitted as gravitational 
waves, these sources are observable to a distance of 
13.2 pc. Conversely, a source at a distance of 15kpc re- 
quires an energy of 1.3 x 10 _5 Mqc 2 to be emitted as 
gravitational waves. This is several orders of magnitude 
greater than the typical energy we might expect from a 
pulsar glitch, making these sources a more attractive tar- 
get for next-generation detectors such as advanced LIGO 
[25l ] and GEO-HF (26j. If, for example, we consider ad- 
vanced LIGO with a factor ~ 10 improvement in sen- 
sitivity over LIGO, the energy required to produce the 
same signal-to-noise ratio drops to ~ 1.0 x 10~ 7 Mqc 2 . 
While this is still well below our fiducial pulsar glitch 
gravitational wave energy of 10 Mqct, it is compa- 
rable to the gravitational wave energy expected to be 
emitted following the axisymmetric collapse of the core 
of a massive star, some fraction of which will be chan- 
nelled into the oscillatory modes of the proto-neutron 
star (27l I2II ]. Again, if we assume the fiducial gravita- 
tional wave energy of 10 _11 AfQC 2 , advanced LIGO will 
be sensitive to neutron star ring-downs at a distance of 
148 pc. 

We find the simple comparison between white noise 
and ring-downs is insufficient to discriminate between 
different types of signal and will be fooled by any tran- 
sient departure from Gaussianity in the noise. Robust- 
ness is improved by including a toy catalogue of glitch 
models (a sine-Gaussian). The odds ratio may then be 
extended (i.e., O123) to mitigate the effects of these un- 
wanted signals at a relatively small sensitivity cost: a 
ring-down signal-to-noise ratio of p = 8.0 is now re- 
quired for the 50% ring-down detection efficiency and 
1% false alarm probability due to Gaussian white noise, 
so that the total gravitational wave energy required for 
50% detection probability of a source at 15kpc is now 
3.3 x 10 _5 Mqc 2 and the observable range is now 8.2 pc. 
With advanced LIGO sensitivity, these estimates improve 
to 2.7 x 10 _8 Mqc 2 and 91.6 pc for the energy in gravita- 
tional waves and the observable range, respectively. We 
also find that this extened odds ratio only has a 10% 
chance of falsely identifying a sine-Gaussian as a ring- 
down. This assumes a worst-case scenario in which the 
sine- Gaussian has a signal-to-noise ratio p ~ 4 (where 
the algorithm's response to these signals peaks). This is 
therefore an upper limit so that the false alarm rate from 
sine-Gaussians will generally be much lower. 

Finally, we have made mention of some of the future 
work required to use the methodology in this work in a 
search for gravitational waves. For the single detector 
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case, it will be necessary to acquire an appropriate cata- 
logue of glitch models and associated prior probabilities 
for their parameters. The extension to the multi-detector 
case will require considerable modification to the wave- 
form model to account for light travel time between de- 
tectors and the appropriate antenna response functions 
for a given sky location. For this search, we are at an 
advantage in that we know the source location and event 
time, and hence the factors introduced by the antenna 
response and time delay are known. We note, however, 
that if this was not so, as would be the case for an all-sky 
search or where we had a trigger but no source location, 
it would be necessary to marginalise over the source lo- 
cation which would significantly complicate the evidence 
calculation. 

With a clear idea of the future analysis pipeline, the 
next step is to characterise the algorithm response using 



real interferometer data and all the consequences of non- 
Gaussian noise. 
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